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GAUGE-INVARIANT FROZEN GAUSSIAN APPROXIMATION METHOD FOR 
THE SGHRODINGER EQUATION WITH PERIODIC POTENTIALS 

RICARDO DELGADILLO, JIANFENG LU, AND XU YANG 


Abstract. We develop a gauge-invariant frozen Gaussian approximation (GIFGA) method for the 
linear Schrodinger equation (LSE) with periodic potentials in the semiclassical regime. The method 
generalizes the Herman-Kluk propagator for LSE to the case with periodic media. It provides 
an efficient computational tool based on asymptotic analysis on phase space and Bloch waves to 
capture the high-frequency oscillations of the solution. Gompared to geometric optics and Gaussian 
beam methods, GIFGA works in both scenarios of caustics and beam spreading. Moreover, it is 
invariant with respect to the gauge choice of the Bloch eigenfunctions, and thus avoids the numerical 
difficulty of computing gauge-dependent Berry phase. We numerically test the method by several 
one-dimensional examples, in particular, the first order convergence is validated, which agrees with 
our companion analysis paper [Delgadillo, Lu and Yang, arXiv:1504.08051 . 


1. Introduction 


The focus of this work is to develop efficient numerical methods for solving the following semiclas¬ 
sical Schrodinger equation whose potential term consists of a (highly oscillatory) microscopic periodic 
potential and a macroscopic smooth potential, 


( 1 . 1 ) 


• A /E 

. 6 ^ = 


Vr (j) + 


X G 


Here is the wave function and e 1 is an effective Planck constant. The equation (1.11 can 

be viewed as a model for electron dynamics in crystal under the one-particle approximation. The 
periodic lattice potential Ur is generated by the ionic cores and electrons in the crystal, and hence 
periodic with respect to the lattice L with unit cell P := [—7r,7r)‘^. In , U is a, smooth external 
macroscopic potential, which counts for e.g., external electric field. 

Direct numerical simulation of ([Oj is prohibitively expensive due to the small parameter e in the 
semiclassical regime. In order to accurately capture the small scale features caused by Up, a mesh 
size of order o(e) is usually required in time and space, e.g., in the standard time-splitting spectral 
method [^. If only physical observables {e.g., density, flux and energy) are needed, one can relax 
the time step requirement to 0{1) with a coarser mesh size of 0(e) using the Bloch decomposition 
based time-splitting spectral method as proposed in [8 TO . However, computation of the solution '0® 
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to (1.1) is still very expensive for e ^ 1, especially in high dimensions. For this reason, alternative 


approaches based on asymptotic analysis have been developed, among which, the geometric optics 
(GO) approach is based on the WKB ansatz under the adiabatic approximation, 

^iS{t,x)/e 


x) = a{t, x)Un VxS, 


Here Un{^, x) is the Bloch eigenfunction normalized for each ^ G F* := [0,1) 
(1.2) j |u„(|,a;)pda; = 1, 

which corresponds to the n-th energy band £'n(4) (see e.g., i): 


d . 


(1.3) 

with the Bloch Hamiltonian 

(1.4) 


1 




Vt{x), 

and periodic boundary conditions on F. 

Then GO solves S{t,x) as the solution to an eikonal equation and p{t,x) = \a{t,x)\‘^ given by a 
transport equation: 

(1.5) dtS + En{V^S) + U{x) =0, 

( 1 . 6 ) dtp + ■ {pV^E„iV^S)) = 0 . 

While this method is e-independent, it breaks down at caustics where the Hamilton-Jacobi equation 


(1.5) develops singularities. 


The Gaussian beam method (GBM) was proposed in to overcome this drawback at caustics, 
with some recent developments [4| [Il||l4{[^[^ , which in particular extends the method to periodic 
media. GBM is based on the single beam solution, which has a similar form as the WKB ansatz 


^/>®(t, a:) = a{t,y)u, 


(V.^, ^ 




S{t,x,y)/e 


The difference lies in that GBM uses a complex phase function. 


(1.7) 


S{t,x,y) = S{t,y) +p{t,y) ■ {x - y) + ]^{x - y) ■ M{t,y){x - y), 


where S' G K, p G M G The imaginary part of M is chosen to be positive definite so 

that the solution decays exponentially away from x = y as a, Gaussian, where y is called the beam 
center. If the initial wave is not in a form of single beam, one can approximate it by a superposition 
of Gaussian beams. The validity of this construction at caustics was analyzed in [^. 

The accuracy of GBM relies on the truncation error of the Taylor expansion of S around the beam 
center y up to the quadratic term, and thus it loses accuracy when the width of the beam becomes 


large, i.e., when the imaginary part of M{t,y) in (1.7) becomes small so that the Gaussian function is 
no longer localized. This happens for example when the solution of the Schrodinger equation spreads 
(the opposite situation of forming caustics). This is a severe problem in general, as shown in 16p9 21 


One can overcome the problem of spreading of beams by doing reinitialization once in a while, see 


20,21 , however, this increases the computational complexity especially when beams spread quickly. 


In the setting of semiclassical Schrodinger equations with periodic potential, another challenge 
for asymptotics methods, not emphasized enough in the literature though, comes from the gauge 
freedom in (1.3). That is, for any Bloch eigenfunction Uni^,x), m„( 4, also solves (1.3) for any 
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arbitrary phase function In particular, when one solves the Bloch waves numerically from (1.3) 

for different it is very difficult, if not impossible, to make sure that the phase depends smoothly on 
The arbitrariness creates difficulty when one needs to get the eigenfunctions off numerical grids by 
interpolation, e.g., in the Gaussian beam method |^. 

In this paper, we develop a gauge-invariant frozen Gaussian approximation (GIFGA) method for the 
Schrodinger equation with periodic potentials. This method generalizes the Herman-Kluk propagator 
by including Bloch waves in the integral representation. It provides an efficient computational 
tool based on asymptotic analysis on phase plane, with a first order accuracy established in our 
companion analysis paper [^. It inherits the merits of the frozen Gaussian approximation studied 
in 16 18 , which works in both scenarios of caustics and beam spreading. The formulation is also 


invariant with respect to the gauge choice of the Bloch eigenfunctions. In particular, we avoid the 
numerical computation of the Berry phase, which causes difficulty since it depends on the derivatives 
of Bloch eigenfunctions with respect to crystal momentum, and is hence not always well-defined if an 
arbitrary gauge choice was made. This is achieved by using a trick inspired by the work of Vanderbilt 


and King-Smith 15 in the context of modern theory of polarization. The details will be explained in 


Section 2.3 see in particular, (2.18|-(2.22). 


The rest of the paper is organized as follows: In Section we will introduce the GIFGA method. 
In Section we briefly describe how to numerically compute Bloch eigenvalues and eigenfunctions. 
We also describe how to numerically implement the GIFGA method described in Section Section]^ 
presents numerical evidence supporting the initial decomposition described in Section along with 
examples conhrming our analytical results in . The last two examples in Section provides the 
numerical performance of GIFGA. We make some concluding remarks in Section]^ 


2. Formulation of the frozen Gaussian approximation 


This section is devoted to the development of the gauge-invariant frozen Gaussian approximation 
(GIFGA) in periodic media based on Bloch decomposition. We first recall the Bloch decomposition 
for Schrodinger operators with a periodic potential. The Bloch waves will be used to capture the high- 
frequency oscillatory structure of the solution given by GIFGA. After stating the asymptotic solution, 
the formulation of which is gauge-invariant, we recall some analytical results on the convergence of 
GIFGA. 


2.1. The Bloch decomposition. Recall that the potential Vr(x) in (l.I) is smooth and periodic 
with respect to the lattice L with unit cell F = [—tt, The unit cell of the reciprocal lattice, known 
as the first Brillouin zone, is then given by F* = [0,1)^^. 

The eigenvalues of the self-adjoint Bloch Hamiltonian H^, defined in (1.4) on T^(r) are real and 
ordered increasingly (counting multiplicity) as 


(2.1) Ai( 0<A2(0<---<A„(0<--- , nSN 

for each ^ G F*. Furthermore, the eigenfunctions {un(^, *)}^i foi' each ^ G F*, known as the Bloch 
waves, form an orthonormal basis of L^(r) [^. 

We extend u„(^,a;) periodically with respect to x so that it is defined on all of and then the 
Bloch decomposition is given by, V/ G L^(IR‘^), 

1 °° r 

(2-2) f(x) = x)e^^'^(5f)„(0 d|, 
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where the Bloch transform B : L? 

( 2 - 3 ) {BfUi) = ^ 

As an analog to the Parseval’s identity, it holds 


L^(r*)^ is given by 
1 


R'i 


Mn(l,y)e ‘^■^/(y)dy. 


(2.4) 


/ l/(a:)pda; = ^ /" | (S/)„(|) d|. 
-'r* 


We denote the phase space corresponding to one band 

(2.5) n :=R‘^ xT* = \ x G | e T*}. 

Let us define the Berry phase, which will be used later, as 

( 2 . 6 ) AniO = ■)) L^T) ■ 

Here we have used the Dirac bra-ket notation (-I-) in quantum mechanics, i.e., 

{f\9)L^n)= [ fgdy, and {f\g) = [ fgdy, 


where / is the complex conjugate of /. Note that the eigenvalue equation (1.31 and its normalization 
only define ’) up to a unit complex number, in particular, for any function (j) periodic in L* 

(2.7) u„(|, x) = x) 

also provides a set of Bloch waves. This is known as the gauge freedom for Bloch waves. It is known 
that (see e.g., [^) we can choose (j) such that u„(^,a;) is smooth in and then the definition (2.6) 
makes sense. However, different gauge choice might give different values of A„(4), and it is also 
difhcult in numerical diagonalization of the Bloch waves to make sure that the phase dependence is 
smooth. We will come back to this delicacy in the development of numerical algorithms. Note that 
from the normalization condition (1.2), A„(^) is always a real number. 

2.2. Formulation. We denote G® ^ the semiclassical Gaussian function localized in the phase space 
at (q,p): 

( 2 . 8 ) 


Gq^p{x) = exp(-|a:-q|V(2e) + ip- {x-q)/e). 


The frozen Gaussian approximation (FGA) solution to (1.1) with the initial condition ipQ is ap¬ 
proximated by 1^, 


(2.9) iPI.Q^{t,x) = (27r£)3rf/2 X] j^o-n{t,q,p)uu{Pn,x/e)Gl^^pJx)e 


iSn{t,q,p)/e 


X (G'q,pW„(p,-/e)|'!/'o) dqdp. 


The right hand side of ( |2.9[ ) sums over all the Bloch bands. For each n, (Q„(t, q, p), P„(t, q, p)) solves 

rmiltonian /i„ 

= V En{P „), 


the equation of motion given by the classical Hamiltonian /i„(q,p) = En{p) -I- U{q): 

dQ 


( 2 . 10 ) 



= -VUiQJ, 


with the initial conditions Q„(0, q,p) = q and Pn{0,q,p) = P- For simplicity, we shall omit the 
subscripts of gradient whenever it does not cause any confusion. 


4 











In (2.9), Sn{t,q,p) is the action associated with the Hamiltonian dynamics (2.10), given by the 
evolution equation 


( 2 . 11 ) 


^ = P„ • Vp/z„(g„, P„) - P„), 


with the initial condition ^^(O, q,p) = 0. The function a„(t, q,p) gives the amplitude of the Gaussian 
function at time t. With the short hand notations 


( 2 . 12 ) 


dz=dq- idp 


Zn — dz {Qn + i-Pn) j 


the evolution equation for a„ is given by 

(2.13) ^ = -ia„^„(P„)-VP(Q„) + ^a„tr(9,P„V2p„(P„)Z-i)-^a„tr(a,g„V2c/„(g„)Z-i) 

with initial condition a„(0, q,p) = 2"^/^ for each {q,p). Recall that yln(4) is the Berry phase of the 
n-th Bloch band given in (|2.6|). 


2.3. Gauge-Invariant Integrator. The gauge freedom of the eigenfunction ri„(^,x) of (1.4) causes 
problems for numerical computation. In particular, different choice of gauge may lead to different 
numerical results for the Berry phase term .4„(^) = (it„(^,a;)|iV^it„(^,a;)), and hence different V'fga 
which is artificial. It is desirable hence to design an algorithm that is manifestly independent of 
the gauge. The key is to avoid direct computation of the the Berry phase and so to avoid the the 
computation of the momentum-gradient of 

First, we separate the dependence of a„ on An in the evolution equation (2.13). For this, we define 
S:^ the phase contribution due to the Berry phase term 

(2.14) 


Let 


.4„(P„) • V[/(g„)ds. 

bn{t,q,p) = a„(t,q,p)exp(iS';(^(t,q,p)), 


then it solves 

(2.15) ^ = Untr(dzPn^^EniPn)Z-^) - tr [dzQn^^U{Q JZ-^) , 

with initial condition bn{0,q,p) = 2“^/^. The evolution equation (2.15) for bn is manifestly gauge- 
invariant, as all terms are independent of the gauge choice. Using the amplitude function the 
frozen Gaussian approximation can be rewritten as 


(2.16) ipl,Qp^{t,x) = 


1 


(27re)3'^/2 7r* 


E 


bn{t, q, p)Un {Pn, X/s) ^x)ASn{i,<i,v)/ 

X (Gq,pMn(p, •/e)IV'o)dqdp. 


The gauge-dependent term in (2.16) thus reads 

(2.17) Un{Pn, £c/e)e“‘'®- y/e). 


Our goal is hence to design a gauge-invariant time integrator for ( |2.14 | such that the term (2.17) 

ge. Observe that, by the Hamiltoniai 

S^it,q,p)=- f A{Pn) ■ dPnis). 

Jo 


becomes independent of the gauge. Observe that, by the Hamiltonian flow (2.10), 

c* 


(2.18) 














Let 0 = to < ti < ■ ■ ■ < Ik = t he a. time discretization, we have 

(2.19) exp(-i5';f) =exp^i J A{Pn) ■ dP„{s)^ = -^(Pn) ■ dP„(s)^. 

To proceed, let us first work in a gauge where m„(^, •) is smooth in ^ € T*. Note that since our final 
formula is gauge-independent, the choice of the gauge here is only for the derivation. Using the Taylor 
approximation, we obtain 

ptk 

i / A{Pn) ■ dP„(s) = -idm{(M„(P„(4_i),-)|VM„(P„(4_i),-)) ■ APfc_„}-HC>(APfc_„)2 


( 2 . 20 ) 


'ik-1 


= iam{l - (■U„(P„(4_i), ■)\Un{Pnitk), •))} + 0{APk,nf 
= i'3m{\n{uniPn{tk), ■)\Un{Pnitk-l), •))} + 0{APk,nf, 

where AP^ „ = P„{tk) — Pn{tk-i)- The first approximation was obtained by using a left Riemann 
sum. The next approximation is the forward difference approximation for the derivative. The last 
approximation is the Taylor series for In 2 : around z = 1. Therefore, taking exponential, we get 


( 2 . 21 ) 


expi 1 


^(p„) • dp„(s)) = + oiAPk,nr. 


tk-l 


I {'^n(P n{tk')i ‘)I^?^(T^ n{tk—l)-i ‘ 


Substituting the last equation in the right hand side of (2.19) gives an approximation to exp(—i^^f) 
with and error 0{APn) with AP„ = max |APfc^„|. This then gives the approximation to (2. as 


(2.22) Un{Pn,xle)e y/e) = P„(<, q,p, a:, y)-f C>(AP„) 


K 


•— I(Pn (^ic) 5 ^/^)} 


/c=l 


{'^n (-^n (^/c) 5 *) 

\'^n{Pni'^k—l') ■! ')) 

1 {Pn (^fc) 5 ') 

\'Un{P n(^/c— 1)5 ')) 1 


(u„(P„(fo),y/e)| -f e>(AP„). 


The right hand side of (2.22) is manifestly gauge-invariant, as the phase term in |Mn(P„(tfe), •)) will 
cancel with that of (u„(P„(tfc), •)!, for k = 0,... ,K. 

Therefore, in summary, we arrive at a gauge-invariant reformulation of V'fga 

(2.23) ^FGA(i:a:) « 715-4^ ^ / f bn{t,q,p)Fn{t,q,p,x,y)G‘^Q p {x) 

X FSAt,<i.P)/e (G- pIV^o) dqdy, 

where is given by (2.22), and the evolution of {Q^,Pn) follows the Hamiltonian dynamics 

dQr, 


(2.24) 


dt 

dl\ 

dt 


— '^En(P „), 

= -vc/(gj. 


with initial condition Qn{h,q,p) = q and P„(0, g,p) = p. 
The action Sn solves 

dSr.. 


(2.25) 


dt 


= Pn ■ Vp/l„(Q„, P„) - /l„(g„, P„), 


with initial condition SniO, q,p) = 0, and the amplitude bn follows the evolution 
(2.26) ^ = ^bntr {d,PnV^EniPn)Z-^) " ^^ntr [/„(Q„)Z-i) , 













with initial condition bn{0,q,p) = 2'^/^. 

2.4. Analytical Results. To make the presentation self-contained, we briefly recall here the analyt¬ 
ical results proved in for the frozen Gaussian approximation to (1.11. The proofs of these results 
and more details can be found in [^. 

First we recall that the FGA ansatz recover the initial condition at time t = 0, '0 fga(Oj = V'o(a^)- 
This follows from the Bloch decomposition (2.21. 

Let us recall a few notions from to state the convergence results for the frozen Gaussian approx¬ 
imation. We define the windowed Bloch transform (W/)„(q,p) : as 


(2.27) 
where 

(2.28) 


{yvf)n{q,p) = 


2dli 

(27r)3d/4 


J'n{P,-)Gq,p\f) = 


2^/4 


(27r)3d/4 


Unip, x)Gq^p{x)f{x) dx 


G. 


q,p(x) : — exp I 


The adjoint operator W* : L^(fl)” — 

2^/4 

(2.29) ' 


\x-q\ 

2 

is then 


+ ip-{x-q) 


20/4 r r 

(W*g)ix) = j^-^^Y^Jj^Unip,x)Gq^p{x)gn{q,p)dqdp. 

The windowed Bloch transform and its adjoint have the following important property. 
Proposition ([^ Proposition 2.2]). The windowed Bloch transform and its adjoint satisfies 
(2.30) W*W = Idi2(R.). 

Remark. Similar to the windowed Fourier transform, the representation given by the windowed Bloch 
transform is redundant, so that WW* ld 2 , 2 (Q)N. The normalization constant in the definition of W 
is also due to this redundancy. 


The previous proposition motivates us to consider the contribution of each band to the reconstruc¬ 
tion formulae (2.30). This gives to the operator : L^(IR'^) —L^(IR'^) for each n € N 

ndji 

(2.31) ' 


nd/A P P 

(n^f){x) = ^2Tr)3d/4 JJ^Un{p,x)Gq.p{x)(Wf)n{q,p)dqdp. 


(27r)3'^/4 

It follows from ( 2.30| ) that 

Correspondingly, the semiclassical windowed Bloch transform : Lfi 


is defined as 


(2.32) {W^fUq,p) = 


2d/4 


(27re)3^/^ 

Similarly we also have the operator : L? 

2dli 


(P,-/e)Gq,p|/) = 


2d/4 


u„(p, x/e)Gg p{x)f{x) dx. 


' (27r£)3^/4 V 

for each n G N with semiclassical scaling 


(2.33) (n.^’’^f){y) = ^ 2 ^^^ 3 d /4 y/£)G%,^iy)(W^f)nix, |) dx d^. 

It follows from ( 2.30| ) and a change of variable that — IdL 2 (Rd). 

For the long time existence of the Hamiltonian flow (2.10), we will assume that the external potential 
U{x) is subquadratic, such that ||9“17(a;)||ic=o is finite for all multi-index jaj >2. As a result, since 
the domain of ^ is bounded, the Hamiltonian is also subquadratic. V’fga provides an approximate 
solution to equation o to first order accuracy as stated in the two theorems below, rephrased from 
our previous work [^. 
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Theorem (|^ Theorem 3.1]). Assume that the n-th Bloch band En{^) does not intersect any other 
Bloch bands for all ^ G F*; and moreover, the Hamiltonian hn{x,^) is subquadratic. Let be the 
propagator of the time-dependent Schrodinger equation (1.1). Then for any given T, 0 < t < T and 
sufficiently small, e < Eq, 

(2.34) sup [ bnit,q,p)ur,{Pn,x/e)Gl^ p {x)x 

0<t<T [fX-KsY^I^ Jq 


X e 


iSrL{t,q,p)/e—iS:^{t,q,p) i/^e 


(Gq,pU„(p, •/£)|V’o)dqdp < Ct,™ e ||^/’§||^ 2 - 


Theorem ([^ Theorem 3.2]). Assume that the first N Bloch bands En{^), n = I,-- - ,N do not 
intersect and are separated from the other bands for all ^ € F*; and assume that the Hamiltonian 
hn{x,^) is subquadratic. Let be the propagator of the time-dependent Schrodinger equation (1.1). 
Then for any given T, 0 < t < T and sufficiently small e, we have 


(2.35) sup 

0<t<T 


1 ^ f 

/ bn{t,q,p)un(Pn,x/e)Gl^ .p„(a:)e 


X (Gq,p'u„(p, •/£)l^/io)dqdp 


\Sn{t,q,p)/e—iS':^{t,q,p) . 


N 


L2 


< Ct,n£\\YI\\l2 + IIV'o - 


n=l 


These approximation results show the first order asymptotic accuracy of FGA, which will be nu¬ 
merically validated in Section 


3. Numerical implementation 


We will now describe the numerical implementation of the gauge-invariant frozen Gaussian ap¬ 
proximation (GIFGA) method. We will restrict ourselves to one spatial dimension in this paper. For 


one thing, the computation of true solutions to (1.1) with high accuracy is extremely time-consuming 


in high dimensions, and thus it is difficult for us to confirm numerically the asymptotic convergence 
order with the pollution of non-negligible numerical errors. For another thing, band-crossing is quite 
common in high dimensional cases (e.g., in honeycomb lattice), which requires more techniques than 
the scope of this paper, and we will leave the numerical study of high dimensional examples as future 
work. The calculation of the Bloch eigenvalues and eigenfunctions is discussed in Section 3.1 In 


Section [372] We describe the numerical algorithms of GIFGA based on the Bloch bands. We will also 
discuss the mesh sizes required for accurate computation. 

3.1. Numerical computation of Bloch bands. We show how to compute numerically the eigen- 

1 d=l. ] 

Un{£.,V) = 


values and eigenfunctions of (1.3) in d = 1. Define the Fourier transform of Un{f,,x) as 

1 


(3.1) 


27T 




'dx. 


Taking the Fourier transform of (1.3) one obtains 

(v + 0^ 


(3.2) 


Un{f,v) + Vr{r]) * ufif„ri) = Enif)un{^,r]), 


where stands for the operation of convolution. 



















Truncating the Fourier grid to{—A,-- - ,A—l}cZ gives 




f Un{i,-A) ^ 


( Unif,-A) ^ 

(3.3) 

H^{A) 

(C; 1 A) 

= En{0 

(^; 1 A) 



Vu„(^,A- l)y 


\Uni^,A- 1)J 


where i^{(A) is the 2A x 2A matrix given by 

/(-A + C)" - 


(3.4) H^{A) = 


Vril) 


VriO) yr(-l) 

(A +1 + ^)2 


■4A(0) 


^ yr(2A-l) V{2A-2) 


yr(l-2A) ^ 
Vr{2-2A) 


(A 


1 + 0^ 
2 


+ VriO)J 


After diagonalizing the matrix, the eigenfunction in the physical domain is then obtained via inverse 
Fourier transform 


(3.5) 


A-1 

V=-A 


Example 3.1. In this example, we compute Bloch eigenvalues and eigenfunctions with potential 
hr (a;) = exp(—25a:^). The extension of Vr{x) periodically with respect to F is not analytic on the 
boundary ofT. However, this lack of smoothness presents a negligible problem numerically as Vr{x) 
decays rapidly. Figure^shows the energy eigenvalues En{f) for ^ € [0,1). The plot shows the first 8 
bands where the bottom curve corresponds to n = 1 (lowest band) and the top curve represents n = 8 
(highest band). Figure^ shows the modules of the corresponding Bloch eigenfunctions for the first 4 
bands. Notice that while these surfaces are continuous and periodic, the next two figures (^ondgj) of 
the real and imaginary parts of the Bloch eigenfunctions are not. This is due to the arbitrary gauge 
freedom in the diagonalization. 


Remark. 1. In the numerical computation of E{fO, the corresponding eigenfunctions and their deriva¬ 
tives near the points ^ = 0 and f — 0.5 (and ^ = 1 by periodicity) is tricky, since the Bloch bands are 
close to each other near these points (see Figure [^. For this reason, our grid for the ^ variable will 
not contain these points. In other words, we shift the grids in the first Brillouin zone to avoid these 
high symmetry points. 

2. One can apply the same technique to derive an algorithm for computing Bloch eigenvalues 
and eigenfunctions in higher dimensions. The main issue with this algorithm is that the numerical 
cost increases drastically for d > 1. In the case where the periodic potential has the form Vt{x) = 
with Vj{xj+2TT) = V{xj), computation of Bloch bands can be treated for each coordinate 
Xj separately. For some common potentials, data for the energy eigenvalues has already been produced 
(see remark 2.1 in [^). 

3.2. Algorithms for gauge invariant frozen Gaussian approximation. We assume that the 
initial data ijjo{x) has compact support or that it decays sufficiently fast as |a:| —> oo, and hence, we 
only need to use a finite number of mesh points in physical space. 
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Bloch Eigenvalues for the first eight bands 



Figure 1. Energy eigenvalues for the one-dimensional lattice potential V (x) 


exp (—25x^) 


For a mesh size Sx and starting point x^ G K, the grid is specified as 
(3.6) x"^ = x'^ + (m — l)Sx, 


for m = 1, • • • , Nx, where is the number of the spatial grid in one dimension. 
We present the algorithm in five steps below. 


Step 1. Compute the Bloch eigenvalues En{^) and eigenfunctions of (1.3), according to the 

algorithm described in Section [XT] 


Remark. For our one dimensional examples in Section]^ we choose a mesh for (^,a;) such that 5^ = 
(1 — 2p)/199 with = —1/2 -I- p and = 200; and 5x = 27r/804 with x^ = —tt and = 805 for 

some 0 < p ^ 1. p was included to avoid putting mesh points at high symmetry points in the first 
Brillouin zone. This number of grid points is enough to ensure that the eigenvalues and eigenfunctions 
are computed with sufficient accuracy for our numerical tests. 


Step 2. Compute (Q„(C ( 7 ,p), P„(C g,p), g,p), 5„(t, g,p)) in ( |2.24| ), ( |2.25[ ), and ( |2.26D . 

To integrate the ODEs for (Qn, S'„, 6^), we use a symplectic fourth order Runge-Kutta method. 
Coefficients for the Butcher tableau can be found in 23 . We will choose a mesh for {q^p) G fl and 
(Qn, Pn) takes initial value at the grid points. That is, 


(3.7) Qn{Q,q,p) = =q° + I5q 

(3.8) Pn{0,q,p) = p-^ =p° + JSp 
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Figure 2. Module of eigenfunctions for the one-dimensional lattice potential V(x) = 
exp (—25x^). We display absolute value of the first 4 lowest energy eigenfunctions. 

where / € I,-- - ,Nj and J G I,-- - ,JVj- Notice that to represent the initial condition ^p(j^(0,a;) 
one only needs the mesh points near x. To be more precise, as the standard deviation of the 
semiclassical Gaussians in ( |2.8| ) is ^/e so one only needs the mesh points contributing significantly 
to 'i/'fga(0j x) satisfy \x — q^\ < This implies that one can put a finite number of mesh points 

for g-coordinate and not on all of M. The mesh size for q^ and p'’ is chosen to be 0{y/£), which 
resolves the oscillation of the initial condition. 

Step 3. Compute the windowed Bloch transformation of the initial condition (m„(p, ■je)G'^^ pIV^o)- For 
the sake of convenience, denote this term by w^[q,p). Let 

(3.9) + {K - l)5y 
be a discrete mesh of y where K = 1, - ■ ■ , Ny. Then, 

Ny 

(3.10) <(<z^p^)« g Gly^yyiy^Mp\y^/e)My^)r9{\y^-q^\)Sy, 

K=1 

with rg a cut-off function such that rg = I'm the ball of radius 0 > 0 centered at the origin and rg = 0 
outside the ball. 

The mesh y^ should approximately cover the support of the initial condition tfjo(y). As can be seen 
by the form of the size of Ny will depend on e. The mesh should be fine enough to accurately 
capture Un{p,y /s)Gy p(y)'ipo{y) for all bands n. 
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Figure 3. Real part of the eigenfunctions for the one-dimensional lattice potential 
V{x) = exp (—25a:^). We display the real parts for the first 4 lowest energy eigen¬ 
functions. We use 100 data points for the | variable. 


Remark. One can reduce the computation time of by incorporating the periodicity of 

Un{^,x) with respect to x. As can be seen by Figure]^ Un[i,x) tends to become more oscillatory as 
n increases. Thus, the mesh of should be adapted so that it depends on n. 


Step 4. Denote the product term in (2.22) by 


(3.11) 


K 

Fn{t,q,p) ■■= 

/c=l 


{u{Pn{tk), ■), u{Pn{tk-l), •)) 
{u{Pn{tk), ■), u{Pn{tk-l), ■)) I ’ 


and note that 

Fnit,q,p,X,y) = \un{PnitK),x/e))Fn{t,q,p){un{Pn{to),y/e)\. 

At this point we now have the required data to compute Fn. Discretize F^ using the same mesh 
from the previous steps to obtain q^,p'^). Here, = 0 <ti < t 2 < ■ ■ ■ <tK = t the temporal 

mesh used in Step 2, with 

t 

^ K’ 
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tj = jSt, St 


and j = 1, ■ ■ ■ , K. 


























Figure 4. Imaginary part of the eigenfunctions for the one-dimensional lattice po¬ 
tential V{x) = exp (—25a:^). We display the imaginary parts for the first 4 lowest 
energy eigenfunctions. We use 100 data points for the ^ variable. 


Step 5. Reconstruct the solution using (2.231 


(3.12) 


N 

~ ,p'^)Un{Pn{t,q^ 


n=l I J 


X Fr,{t,q^,p-’)il}l{q\p-’)rg {\x^ - QI^'’\)^5q5p, 


where Q„ and P„ are evaluated at (t, q^,p^), and rg is a cutoff function as described in Step 3 and N 
is the maximum number of Bloch bands used. 


Accuracy. The theorems in Section |2.4| and (2.22) imply the above algorithm has a total accuracy 
O (e -|- -|- max„ where 0(6t'^/e) comes from the approximation 

to the phase functions in ( |2.23 ) and 0(max„AP„) « 0{5t) is due to the approximation (2.22). 

' 0 Q||r ,2 is the initial decomposition error, which in general decays with the number 
of bands as indicated in, e.g.. Examples 4.1 and 4.2 and in |^. 
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Table 1 
of £ and 


£ = 1/64 

Error \\po - '/'fgaIU^^ 

iV = 1 

0.13260 

N = 2 

0.11328 

N = 4 

0.033126 

N = 8 

7.2587e-05 

£ = 1/128 

Error \\po - ^gaIU^ 

iV = 1 

0.15361 

N = 2 

0.096905 

N = 4 

0.031652 

II 

oo 

7.0574e-05 

£ = 1/256 

Error \\po - PfqaWl^ 

N = 1 

0.14165 

N = 2 

0.1063 

N = 4 

0.032405 

II 

oo 

6.9192e-05 

£ = 1/512 

Error \\po - 

iV = 1 

0.15885 

N = 2 

0.09276 

N = 4 

0.031263 

N = 8 

6.8701e-05 


. error of ipoix) — V'fga( 0; for Example 
sum over N Bloch bands in V'fga- 


4.1 


We display various values 


4. Numerical examples 


In this section, we show the numerical performance of gauge invariant frozen Gaussian approxi¬ 
mation (GIFGA) by several one dimensional examples, which also confirm the first order asymptotic 
convergency analyzed in [^. 


4.1. Initial decomposition. In the first two examples, we test the initial decomposition of GIFGA 
described in Section]^ We compute '0 fqa at t = 0 via equation (2.9). As we cannot numerically 
sum to infinity, we choose to use at most 8 bands in all of our examples. Expressed differently, the 
solution will be concentrated on the first 8 bands. Because of the need for 0{^/e) mesh size for both 
coordinates {q^,p'^) of phase space, we choose approximately 2/^/e number of grid points for each 
unit interval. 


Example 4.1. In this example, we check the initial decomposition by choosing = A(a;) exp(iS'(a:)/£) 
with A{x) = exp (—50a;^) cos((a; — 0.5)/£) and S{x) = 0.3(a; — 0.5) -I- 0.1sin(a; — 0.5), and the lattice 
potential Vr = cos(a;). We record the data in Ta 6 Ze[^ 


Example 4.2. In this example, we check the initial decomposition by choosing po = A(a;) exp(iS'(a:)/£) 
with A{x) = exp (—50a;^) and S{x) = 0.3 -I- 0.1sin(a; — 0.5), and the lattice potential to he Vr = 
exp(—25a:^). We record the data in Table^ 
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£ = 1/64 

Error Hi/o - V'fga 

iV = 1 

0.035736 

N = 2 

0.02463 

N = 4 

0.0075756 

N = 8 

0.0018796 

£ = 1/128 

Error \\ipo - V’fga! 

iV = 1 

0.031445 

N = 2 

0.024814 

N = 4 

0.007579 

II 

oo 

0.0018579 

£ = 1/256 

Error \\ipo - V’fgaII i’’ 

N = 1 

0.030633 

N = 2 

0.024967 

N = 4 

0.0076045 

II 

oo 

0.0018698 

£ = 1/512 

Error ||V’o - V’fga! li’’ 

iV = 1 

0.030375 

N = 2 

0.025078 

N = 4 

0.0076103 

N = 8 

0.0018769 


Table 2. error of '0o(a^) — V'fga(0; for Example 
of £ and sum over N Bloch bands in V'fga- 


4.2 


We display various values 


Tables and show that EGA indeed matches the initial condition more closely as N increases. 
Eurthermore, we have essentially the same error for each e. This provides numerical verification 
of the independence of e of the initial decomposition. 


Remark. Let us note that from equation (1.4) the convergence rate should depend on the form of the 
lattice potential Vr(a;)- Also, by equation (2.2), the convergence rate also depends on the form of the 
initial condition. We see from Examples 4.1 and 4.2 that the cosine lattice potential seem to produce 
faster convergence with respect to the number of bands used. Different initial conditions may also 
converge faster as N increases. Example |4.4| uses an initial condition projected onto the first band. 
Choosing such initial condition has the advantage of needing only to compute V'fga band. 

By examining the errors or the convergence rates, one could determine the minimum number 
of bands to sum over to achieve required accuracy. In Example |4.1[ it shows that upon summing over 
iV = 4 bands, the initial decomposition starts to resemble the initial condition. 


4.2. Verification of the convergence rate of EGA. Eirst, we choose to test the convergence rate 
of (2.9) with external potential U(x) = 0 in Examples |4.3|and|4.4[ With this choice of potential, there 


is no need for a gauge-invariant algorithm. One can optimize the algorithm described in Section 3.2 


by setting F{t,q,p) = 1 in ( |3.11[ ) in Step 4. Thus, for Examples |4.3| and 

and 


from F{t,q,p) will be absent. Examples 


4.5 


4.6 


4.4 


numerical errors coming 


have nonzero external potential so there will 
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N=1 



N=2 





Figure 5. Initial decomposition for example 4.1 The real part of and 


V'fga( 0>*) are shown for e = 1/256. The summation in '(/;pQ^(0,x) is over the 
first 4 lowest energy bands. 



Figure 6 . The plot of ||V'o(a;) — '0 |qa(O,cc)||/2 for figure]^ is displayed here. 


be some numerical errors introduced by F{t,q,p). We continue using ^2,1^ mesh points per unit 


interval in q and p and sum up to eight bands (except for Example 4.4). We choose a time step of 


size At = T/150. The exact solution to equation (1.1) will be computed using the Strang splitting 


spectral method [^. For all of our examples, the Strang splitting spectral method did not need a 
mesh finer than Aa: = 1/2^® and At = 1/2^^. 
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Error Wifspec - V’fgaIU^ 

Rate of Convergence 

e = 1/8 

0.09112 


e = 1/16 

0.048907 

0.8977 

e = 1/32 

0.022603 

1.1135 

e = 1/64 

0.010555 

1.0986 


Table 3. error of '0spec(O.35, x) — '(/'FGA(0-35,a::) for various values of e. The 
summation in V'fga over the first 8 lowest energy bands. 




4.3 


Figure 7. Example 
side with the error for e = 1/8. 


plot of real parts of '0 |qa(O. 35, x) and ^spec(0.35, x) along 


Example 4.3. In this example we choose the initial condition to be tpo = A(x) exp(i5'(x)/£) with 
A{x) = exp (—50x^) and S{x) = 0.3 + 0.1 sin(x —0.5). The exact solution is computed using the Strang 
Splitting spectral method. This is done at time T = 0.35. The lattice potential used is Vr(a^) = cos(x). 
We record the data in The convergence order of the data in table^is 1.0366. ITe display 

plots of the solution for e = 1/8,1/16,1/32 and 1/64 in Figures and \10l[ 


In the next example, we will choose initial condition projected onto one Bloch band. With this 
choice of initial condition, there will be no initial error. 

Example 4.4. In this example we will choose an initial condition Il^ff'ipo^x) given by (2.33) with 
tpo{x) = 7l(x) exp(i5'(x)/e) where A{x) = exp(—50x^) and S{x) = 0.3x + 0.1sin(x — 0.5) with lattice 
potential exp(—20x^) and external potential U(x) = 0. We compute the solution at time T = 0.35 
using the Strang Splitting spectral method and GIFGA. The errors are recorded in The 

convergence order is 0.9814- We display plots of the solution for e = 1/64,1/128 and 1/256 in 
Figures [77} [7^ and \13\ 
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X X 


Figure 8. Example 4.3 plot of real parts of '0fqa 


side with the error for e = 1/16. 


(0.35, x) and ^5pec(0.35, x) along 



X X 


Figure 9. Example 4.3 plot of real parts of i/fga 


side with the error for e = 1/32. 


(0.35, x) and ^spec(0.35, x) along 
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4.3 


Figure 10. Example 
side with the error for e: = 1/64. 


plot of real parts of '0pQ^(O.35, a;) and '(/'Spec(0.35, x) along 



Error Wifspec - 

Rate of convergence 

e = 1/64 

0.0269 


e = 1/128 

0.0144 

0.9015 

e = 1/256 

0.0069 

1.0614 


Table 4. error of i/>spec(0.35,x) — tj;pQj^^{0.35,x) for initial condition projected 
onto the first Bloch band. 


Example 4.5. In this example we choose the initial condition to be tpo = A(x) exp(i5'(x)/£) with 
A{x) = exp (—50x^) cos((x — 0.5)/e) and S{x) = 0.3(x — 0.5) + 0.1 sin(x — 0.5). The exact solution is 
computed using the Strang Splitting spectral method. This is done at time T = 0.2. The potential used 


is Vr(x) = exp(—25x^) with external potential U{x) = 




Our results are shown in Table 


convergence order of the data in 0.9488. We display plots of the solution for e = 1/128 

and 1/512 in Figures I 4 . 15, and\l(^ 


The 
1/256 


Example 4.6. In this example we choose the same initial condition as in Example \4.5[ All of the same 
parameters as in Example \4.5\ will also be used. The exact solution is again computed using the Strang 
Splitting spectral method at time T = 0.2. The only difference is that we change the external potential 
to U{x) = cos(x). The convergence order of the data in ra6Ze[^ is 0.8439. We display plots of the 
solution for e = 1/128,1/256 and 1/512 in Figures [T^ [7^ and \ 1 
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Figure 11. Example 4.4 plot of the real part of 4’Speci0.35,x) and 4 ’fga 

b'pec(0.35, a;) — '0|ng.^(O.35, x) for example 


alongside with the error of ipsp 
£ = 1/64. 


4.4 


(0.35, x) 

. We use 



X X 


4.4 


Figure 12. Example 
alongside with the error of ipS: 
£ = 1/128 . 


plot of the real part of '0Spec(O.35, x) and 4 ’'fga 
■pec(0-35, x) —'!/)|.g,^(0.35, x) for example 


4.4 


(0.35, x) 

, We use 
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Figure 13. Example 4.4 plot of the real part of 4’Speci0.35,x) and 4’fga 

b'pec(0.35, a;) —x) for example 


alongside with the error of ipsp 
e = 1/256 . 


4.4 


(0.35, x) 

. We use 



Error \\ipspec - V’fgaHl^ 

Rate of Gonvergence 

e = 1/64 

0.059576 


e = 1/128 

0.038811 

.61826 

e = 1/256 

0.015225 

1.3500 

£=1/512 

0.0082833 

0.8782 


Table 5. error of 4’Spec{0-2,x) — '0J.g,^(O.2, x) for various values of e. The sum¬ 
mation in ijjpQA is over the first 8 lowest energy bands. 



Error \\ipspec - ^I^fgaWl^ 

Rate of Gonvergence 

£ = 1/128 

0.039714 


£ = 1/256 

0.019057 

1.0593 

£ = 1/512 

0.012327 

0.6285 


Table 6. error of ijjSpec{0-2,x) — i/’f.g^(0.2, x) for various values of e. The sum¬ 
mation in 4’fga i® over the first 8 lowest energy bands. 


5. Discussion and Conclusion 

In this paper, we generalize the Herman-Kluk propagator for the linear Schrodinger equation (LSE), 
and develop the gauge-invariant frozen Gaussian approximation method for LSE with periodic po¬ 
tentials in the semiclassical regime. The method is invariant with respect to the gauge choice of the 
Bloch eigenfunctions, and thus avoids the numerical difficulty of computing gauge-dependent Berry 
phase. The numerical examples show that that the frozen Gaussian approximation is indeed a good 
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X X 


4.5 


Figure 14. Example 
along side with the error for e = 1/128. 


plot of the real parts of '0pQ^(O.2, *) and 'ipspec{^-‘^,x) 



- 11 (0.2, CC) - «^FG 4 ( 0 .2 , x) 11 


0.25 


0.2 


0.15 


0.1 



Figure 15. Example 4.5 plot of the real parts of %l)pQj^{Q.2,x) 


along side with the error for e = 1/256. 


and i/spec(0.2,a;) 
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0.8 


0.6 - 


---f^e['^s'pec(0-2, 

a;)] 



- 11 (0.2, x) - 4^fg.4(0.2, x) 11 

- Re\^FGA{^-'^ 

a;)] 


0.25 




0.2 


0.15- 


0.1 - 


0.05- 


-.5 -.25 



Figure 16. Example 4.5 plot of the real parts of '0 |.q^(O. 2, *) and 'ipspeciQ-‘2',x) 
along side with the error for e = 1/512. 



X X 


Figure 17. Plot of real parts of 'tljpQj^{0.2,x) and ipspec{0-‘^,x) along side with the 
error for e = 1/128. 
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Figure 18. Plot of real parts of '0 |^q^(O. 2, x) and ipspec{0-‘^,x) along side with the 
error for e = 1/256. 



X X 


Figure 19. Plot of real parts of 'tljpQj^{0.2,x) and ipspec{0-‘^,x) along side with the 
error for e = 1/512. 
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approximation to the exact solution of the Schrodinger equation 

order of our numerical results confirms the estimate given in 

(5-1) W-fpExactit^x) - '‘pFGAit,x)\\L 2 = 0(e). 

In future, we will study high dimensional examples where band-crossing happens quite common, and 

thus requires more techniques than the scope of this paper. 
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